Socioeconomics of Pedestrian Accidents
Start by loading the various libraries and packages we will be using.
The “Pacman” package will automatically install any packages you do not currently have installed and load the applicable libraries.
if (!require("pacman")) install.packages("pacman")
pacman::p_load(dplyr, spdplyr, sf, sp, spdep, spatialEco, rgdal, geojsonio, geojsonR, leaflet, leaflet.minicharts, leaflet.extras,leafsync, maptools, tmap, tidyverse, ggplot, gclus, rpart)Create our Data Tables
KSI Pedestrian Dataset
This dataset is a subset of the Killed and Seriously Injured (KSI) dataset collected by the Toronto Police Service from 2008-2018. These events include any serious or fatal collisions where a Pedestrian is involved. To learn more about Pedestrians related collisions in Toronto you can follow this link: http://data.torontopolice.on.ca/pages/pedestrians
We will use the FROM_GeoJson command from the “geojsonR” package to download a json file from the provided URL. The geojson_read command from the “geojsonio” and formating form the “sp” package can then be used to read the json and generate a SpatialPointsDataFrame
ksi_ped <- geojsonio::geojson_read("https://opendata.arcgis.com/datasets/3dedc9bff625450990b8d480f397ad3f_0.geojson", what = "sp")
head(ksi_ped)class : SpatialPointsDataFrame
features : 6
extent : -79.45759, -79.36952, 43.6681, 43.71057 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 55
names : Index_, ACCNUM, YEAR, DATE, TIME, Hour, STREET1, STREET2, OFFSET, ROAD_CLASS, District, LATITUDE, LONGITUDE, LOCCOORD, ACCLOC, ...
min values : 5413936, 1062115, 2008, 2008/09/02 04:00:00+00, 917, 9, REDPATH AVE, BROADWAY AVE, , Collector, Toronto and East York, 43.668101, -79.45759, Intersection, At Intersection, ...
max values : 80632109, 6001446484, 2016, 2016/08/16 04:00:00+00, 920, 9, WELLESLEY ST E, PARLIAMENT ST, , Minor Arterial, Toronto and East York, 43.710574, -79.369525, Intersection, At Intersection, ... KSI TTC/Municipal Vehicle Dataset
This dataset is a subset of the Killed and Seriously Injured (KSI) dataset collected by the Toronto Police Service from 2008-2018. These events include any serious or fatal collision involving an operator or passenger of a TTC, Transit Vehicle, streetcar or Municipal Vehicle. To learn more about TTC-Municipal Vehicle related collisions in Toronto you can follow this link: http://data.torontopolice.on.ca/pages/ttc-municipal-vehicle
We will use the FROM_GeoJson command from the “geojsonR” package to download a json file from the provided URL. The geojson_read command from the “geojsonio” and formating form the “sp” package can then be used to read the json and generate a SpatialPointsDataFrame
ksi_ttc <- geojsonio::geojson_read("https://opendata.arcgis.com/datasets/dc4751278e604d65b0886b9765d4b551_0.geojson", what = "sp")
head(ksi_ttc)class : SpatialPointsDataFrame
features : 6
extent : -79.44409, -79.30189, 43.65105, 43.72284 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 55
names : Index_, ACCNUM, YEAR, DATE, TIME, Hour, STREET1, STREET2, OFFSET, ROAD_CLASS, District, LATITUDE, LONGITUDE, LOCCOORD, ACCLOC, ...
min values : 6597284, 1220073, 2011, 2011/01/06 05:00:00+00, 1404, 14, QUEEN ST W, DUFFERIN ST, , Major Arterial, Scarborough, 43.651045, -79.444089, Mid-Block, , ...
max values : 7645393, 1333336, 2012, 2012/12/31 05:00:00+00, 1415, 14, VICTORIA PARK AVE, YORK ST, , Major Arterial, Toronto and East York, 43.722845, -79.30189, Mid-Block, , ... Merge KSI_TTC and KSI_Ped
As you will have noticed from the descritpions both the KSI_PED and KSI_TTC datasets are themselves subsets of the Toronto Police Services’ KSI Dataset. We downloaded them as seperate dataframes to enable faster downloads as they make up a small portion of the full KSI dataset, which would take significantly more time to download and sort. Because they come from the same base dataset and have the same schema we can merge them back into one dataframe to work with. We will merge them based on the “Index_” to ensure there are no duplicates created due to the merging process.
ksi_merged <- merge(ksi_ped,ksi_ttc, by="Index_")Neighbourhood Boundaries - City of Toronto
The City of Toronto also maintains a list that defines the boundaries of all the neighbourhoods in the city. A file containing the spatial data required to map these neighbourhoods can be downloaded form the City of Toronto open data portal. Due to how this file will download, unlike the KSI data, we cannot directly read the GeoJson file but have to download it before the file can be read.
trying URL 'https://ckan0.cf.opendata.inter.prod-toronto.ca/download_resource/a083c865-6d60-4d1d-b6c6-b0c8a85f9c15?format=geojson&projection=4326'
Content type 'application/json' length unknown
downloaded 1.9 MBclass : SpatialPolygonsDataFrame
features : 6
extent : -79.52174, -79.39119, 43.64744, 43.78728 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 15Simply Analytics - Census Data
From Simply Analytics i have downloaded shapefiles containing average income by census track and subdivision. These two shapefiles will usefull in drilling down beyond the neighbourhood. Both shapefiles can be uploaded into a spatial dataframe using the st_read command from the “sf” package we installed earlier then converting them to a SpatialPolygonsDataFrame.
ogrInfo(dsn="datasets/SimplyAnalytics_C1", layer="C1")Source: "C:\Users\gvs_j\Documents\GitHub\CSDA-1050F18S1\jacobgvs_304292\sprint_2\Datasets\SimplyAnalytics_C1", layer: "C1"
Driver: ESRI Shapefile; number of rows: 572
Feature type: wkbPolygon with 2 dimensions
Extent: (-79.6393 43.56034) - (-79.11347 43.85547)
CRS: +proj=longlat +datum=WGS84 +no_defs
LDID: 87
Number of fields: 4 ctr <- readOGR(dsn = "datasets/SimplyAnalytics_C1", layer = "C1")OGR data source with driver: ESRI Shapefile
Source: "C:\Users\gvs_j\Documents\GitHub\CSDA-1050F18S1\jacobgvs_304292\sprint_2\Datasets\SimplyAnalytics_C1", layer: "C1"
with 572 features
It has 4 fieldsnames(ctr@data)[names(ctr@data)=="VALUE0"] <- "Household Total Income After-Tax"
names(ctr@data)[names(ctr@data)=="VALUE1"] <- "Household Aggregate Income"
names(ctr@data)[names(ctr@data)=="VALUE2"] <- "Household Average Income"
head(ctr)class : SpatialPolygonsDataFrame
features : 6
extent : -79.47429, -79.30083, 43.61037, 43.66591 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 4ogrInfo(dsn="datasets/SimplyAnalytics_C2", layer="C2")Source: "C:\Users\gvs_j\Documents\GitHub\CSDA-1050F18S1\jacobgvs_304292\sprint_2\Datasets\SimplyAnalytics_C2", layer: "C2"
Driver: ESRI Shapefile; number of rows: 3702
Feature type: wkbPolygon with 2 dimensions
Extent: (-79.6393 43.56034) - (-79.11347 43.85547)
CRS: +proj=longlat +datum=WGS84 +no_defs
LDID: 87
Number of fields: 5 disem <- readOGR(dsn = "datasets/SimplyAnalytics_C2", layer = "C2")OGR data source with driver: ESRI Shapefile
Source: "C:\Users\gvs_j\Documents\GitHub\CSDA-1050F18S1\jacobgvs_304292\sprint_2\Datasets\SimplyAnalytics_C2", layer: "C2"
with 3702 features
It has 5 fieldsnames(disem@data)[names(disem@data)=="VALUE0"] <- "Household Total Income After-Tax"
names(disem@data)[names(disem@data)=="VALUE1"] <- "Household Aggregate Income"
names(disem@data)[names(disem@data)=="VALUE2"] <- "Household Average Income"
summary(disem)Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -79.63930 -79.11347
y 43.56034 43.85547
Is projected: FALSE
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
spatial_id name Household Total Income After-Tax
35200002: 1 DA0002, Toronto, ON: 1 Min. : 0.0
35200003: 1 DA0003, Toronto, ON: 1 1st Qu.: 159.0
35200004: 1 DA0004, Toronto, ON: 1 Median : 202.5
35200005: 1 DA0005, Toronto, ON: 1 Mean : 300.6
35200006: 1 DA0006, Toronto, ON: 1 3rd Qu.: 310.0
35200007: 1 DA0007, Toronto, ON: 1 Max. :4889.0
(Other) :3696 (Other) :3696
Household Aggregate Income Household Average Income
Min. : 0 Min. : 0
1st Qu.: 14441832 1st Qu.: 72481
Median : 20666158 Median : 90304
Mean : 30880924 Mean : 116745
3rd Qu.: 33046773 3rd Qu.: 120148
Max. :534324916 Max. :2340057
head(disem@data)Additional Datatables
One additional datatable may be needed depending on how we choose to cluster the KSI data. At a high level grouping accidents to a given neighbourhood as defined by the city then using the census data to further cluster by census track and subdivision will be ideal. I may need to pull the Neighbourhood Profiles data set from the city of Toronto to get an average income level by neighbourhood as this does not appear to be avaialble directly from the census datasets.
Reviewing the Datasets
With our datasets downloaded we can start by looking through the data. Becasue we are dealing almost exclusively with spatial data the easiest way to get a handle on the data is by mapping it. To do this I will be using features from the “leaflet” package.
KSI_merged Dataset
Lets satrt by getting an idea of where accidents are occuring. The leaflet package will map the data for us. To make this more readable I have had the system autocluster the accidents. These clusters don’t have any relation to specific neighbourhoods and will need to be adjusted later so that the clusters are in line with our other datasets.
leaflet(ksi_merged) %>%
addTiles() %>%
addMarkers(lng = ksi_merged$LONGITUDE, lat = ksi_merged$LATITUDE, clusterOptions = markerClusterOptions())This map will allow you to zoom in and the clusters will auto adjust as you zoom in and out. These clusters are based on proximity to a central point. Once you get to the lowest zoom levels, clicking on a cluster will map the individual accidents. Details for each accident are not currenlty included in the mapping.
Neighbourhood Boundaries - City of Toronto
Our next dataset contains the boundary lines for various neighbourhoods in Toronto. mapping this will begin to give some dimension to how we intend to cluster accidents going forward.
leaflet(nbh) %>%
addTiles() %>%
addPolygons()Simply Analytics - Census Data
We can similarly plot both the census tract and dissemination area files pulled from Simply Analytics.
leaflet(ctr) %>%
addTiles() %>%
addPolygons()leaflet(disem) %>%
addTiles() %>%
addPolygons()Our next step will be to assing each accident in the KSI data to a dissemination area polygon. But before we can do that we will first have to clean up the various datasets.
Cleaning the Data
Before we begin we need to ensure that our various datasets are compatible.
Using the CRS fuction in the SF package We can check the coordinate reference system being used by our main files and can see that they are all the same.
st_crs(ksi_merged)Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"st_crs(ctr)Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"st_crs(disem)Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"st_crs(nbh)Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"We luckily all the files we have downloaded use the same co-ordinate reference system so we will be able to compare and associate them to one another without having to re-project them into a comon co-ordinate system. This is especially convenient given that the mapping package w are using “leaflet” does not support “CRS” one of the other reference systems.
Cleaning the KSI data
As we review the datasets one thing you may have noticed is that the KSI dataset contains a large number of potenitally duplicate entries.
head(ksi_merged@data, 20)While each accident is assigned a different index id number(Index_), in most cases multiple index numbers are assigned to the same account number(ACCNUM). This is because, for each accident, the driver of the vehicle, the pedestrian, and any other applicable individual related to the loss has their information recorded to the dataset. This is most evident if you look at the involved vehicle type (INVTYPE) columns.
Since this review is focused on pedestrian accidents, we will filter our KSI_Merged file to include only instances where the column INVTYPE is a pedestrian.
We can easily identify all unique items in the INVTYPE column and select those that are representative of pedestrians.
unique(ksi_merged@data$INVTYPE) [1] Driver Pedestrian Vehicle Owner Passenger
[5] Truck Driver Other Wheelchair Pedestrian - Not Hit
[9] Motorcycle Driver Driver - Not Hit Other Property Owner
[13] In-Line Skater Runaway - No Driver Cyclist Witness
[17] Trailer Owner
17 Levels: Cyclist Driver Driver - Not Hit In-Line Skater Motorcycle Driver ... WitnessIf we look at our list of options, there are a few different entries we should be considering to be a pedestrian. For this review we will be filtering out all entries where the INVTYPE is “Pedestrian”, “Pedestrian - Not Hit”, “In-Line Skater”, or “Wheelchair”.
So as not to lose our merged dataset, we will also create a new filtered datatable. We can easily identify and, if required in the future, modify the list of items being filtered into our datatable by setting them as a target
target <- c("Pedestrian", "Pedestrian - Not Hit", "In-Line Skater", "Wheelchair")
ksi_modified <- ksi_merged[ksi_merged@data$INVTYPE %in% target,]There are also a number of columns that will not be needed for our review and can be removed from our modified dataset. These include: OFFSET, VEHTYPE, MANOEUVER, DRIVACT, DRIVCOND, CYCLISTYPE, CYCACT, CYCCOND, CYCLIST, TRSN_CITY_VEH, TRSN_CITY_, coords.x1, coords.x2
target_coll <- c("OFFSET", "VEHTYPE", "MANOEUVER", "DRIVACT", "DRIVCOND", "CYCLISTYPE", "CYCACT", "CYCCOND", "CYCLIST", "TRSN_CITY_VEH", "TRSN_CITY_", "coords.x1", "coords.x2")
ksi_modified <- ksi_modified[,-which(names(ksi_modified@data) %in% c(target_coll))]
head(ksi_modified@data, 20)Cleaning the Toronto Neighbourhood data
Next we will look at the Toronto Neighbourhood dataset.
head(nbh@data, 20)We can see that for the most part all columns in this data set are useful. The columns that will not be of use and can be removed are: PARENT_AREA_ID - this field containe no unique values AREA_LONG_CODE - this is variable is identical to AREA_SHORT_CODE AREA_DESC - this is variable is identical to AREA_NAME X - this field containe no unique values Y - this field containe no unique values
unique(nbh@data$PARENT_AREA_ID)[1] 49885unique(nbh@data$X)[1] <NA>
Levels: unique(nbh@data$Y)[1] <NA>
Levels: targetnbh <- c("PARENT_AREA_ID", "AREA_LONG_CODE", "AREA_DESC", "X", "Y")
nbh_modified <- nbh
nbh_modified <- nbh_modified[,-which(names(nbh_modified@data) %in% c(targetnbh))]
head(nbh_modified@data, 20)Cleaning the Census Tract and Disemination data
Next we will look at the Toronto Neighbourhood dataset.
head(ctr@data, 20)We can see that these have a very minimal number of fields and will not require any adjustmeents prior to proceeding with our review.
Associating the Datasets
We are interested in associating all our spatial polygon dataframes to the KSI data points that represent each pedestrian related accident.
Our first task will be to join our KSI dataset to a Toronto Neighbourhood.
This can be done using the spatialEco package. The spatialEco::point.in.poly function intersects point and polygon feature classes and adds polygon attributes to points. This function will re-name columns with similar names so we will also address re-naming some of the newly added columns.
ksi_coord <- spatialEco::point.in.poly(ksi_modified, nbh_modified)although coordinates are longitude/latitude, st_intersects assumes that they are planar
although coordinates are longitude/latitude, st_intersects assumes that they are planarnames(ksi_coord@data)[names(ksi_coord@data)=="LONGITUDE.x"] <- "LONGITUDE"
names(ksi_coord@data)[names(ksi_coord@data)=="LATITUDE.x"] <- "LATITUDE"
names(ksi_coord@data)[names(ksi_coord@data)=="LONGITUDE.y"] <- "LONGITUDE.nbh"
names(ksi_coord@data)[names(ksi_coord@data)=="LATITUDE.y"] <- "LATITUDE.nbh"
head(ksi_coord@data, 20)Next we will add in the disemination area information. Being the most granular this will allow us to accurately associate an average income to each accident.
ksi_coord <- spatialEco::point.in.poly(ksi_coord, disem)although coordinates are longitude/latitude, st_intersects assumes that they are planar
although coordinates are longitude/latitude, st_intersects assumes that they are planarhead(ksi_coord@data, 20)Now that we have a assigned an average income to each accident we can start working to review our data and learn more abot what story it can tell.
Since we are trying to look at socioeconomiv impact on accidents we will first start by focusing on average income.
Lets start by creating a histogram of the Household Average Income information for our accidents.
qplot(ksi_coord$Household.Average.Income,
geom="histogram",
binwidth = 5,
main = "Histogram Average Income",
xlab = "Average Income",
fill=I("blue"),
col=I("red"))
As we can see the majority of losses are clustered at the bottom end of the income scale. However this histogram can be deceiving. The average income varies so greatly between the lowest and highest values that it is hard to identify what income brakets this cluster actually represents.
We can re-present the same graph with a reduced bin width but it does little to make the data points more meaningful.
qplot(ksi_coord$Household.Average.Income,
geom="histogram",
main = "Histogram Average Income",
xlab = "Average Income",
fill=I("blue"),
col=I("red"))
To allow us to more easily review the data set we will also convert it from a SpatialPointsDataFrame to a regular dataframe. We can then pull some basic information about household income.
ksi_sf <- as.data.frame(ksi_coord@data)
mean(ksi_sf$Household.Average.Income, na.rm = TRUE)[1] 100023mode(ksi_sf$Household.Average.Income)[1] "numeric"range(ksi_sf$Household.Average.Income, na.rm = TRUE)[1] 0 1065780quantile(ksi_sf$Household.Average.Income, na.rm = TRUE) 0% 25% 50% 75% 100%
0.0 64743.0 82846.0 108044.8 1065780.0 quantile(ksi_sf$Household.Average.Income, na.rm = TRUE, probs = c(0.01, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.75, 1.00)) 1% 2.5% 5% 10% 20% 30% 40% 50% 75%
28453.35 36927.38 44130.50 51818.00 61222.00 68192.00 75706.00 82846.00 108044.75
100%
1065780.00 The quantile results are very interesting in that we can see that the bottom 25% of incomes actually varies between $0 and $65,000 and the highest value is well beyond being plausibly compared to the rest
This does give some indication that we will probably want to focus on the bottom 50% of incomes depedning on how the income level corelates to frequency of accidents
ksi_sf <- as.data.frame(ksi_coord@data)
freq <- table(ksi_sf$YEAR, ksi_sf$AREA_NAME, ksi_sf$ROAD_CLASS, ksi_sf$Household.Average.Income)
yearfreq <- ksi_sf %>%
group_by(YEAR) %>%
summarise(freq = n())
yearfreq <- as.data.frame(yearfreq)
AreaLossFreq <- ksi_sf %>%
group_by(AREA_NAME) %>%
summarise(freq = n())
AreaLossFreq <- as.data.frame(AreaLossFreq)
RoadClassFreq <- ksi_sf %>%
group_by(ROAD_CLASS) %>%
summarise(freq = n())
RoadClassFreq <- as.data.frame(RoadClassFreq)
AveIncomeFreq <- ksi_sf %>%
group_by(Household.Average.Income) %>%
summarise(freq = n())
AveIncomeFreq <- as.data.frame(AveIncomeFreq)A standard frequency plot will do much better at indicating the average income of the areas where these pedestrian accidents usually occur but again will not provide any relevant data
plot(AveIncomeFreq)
We will want to use the quantile results to set income bandwiths to apply to each of our accidents. As evident from the hitograms and frequency plot we will also want to look at focusing on the lower quantiles.
To do:
Establish income bandwiths and associate back to main datatable.
Continue with frequency revie and comparison of various variables.
Correct clustering as part of the evaluation. Example below does not help given it cannot be read.
library(pvclust)
d <- dist(ksi_sf, method = "euclidean")
fit <- hclust(d, method="ward.D2")
plot(fit)
groups <- cutree(fit, k=5)
rect.hclust(fit, k=5, border="red")
- Review how to bind markers in leaflet to their polygon(s). Currently unclear if this feature is available in Rstudio or only available using the leaflet code html or leaflet in python
Lcontent <- paste("Year:", ksi_coord@data$YEAR,
"Income:",ksi_coord@data$Household.Average.Income,
"Location:",
ksi_coord@data$STREET1)
leaflet(ksi_coord) %>%
addTiles() %>%
addCircleMarkers(
stroke = FALSE, fillOpacity = 0.5,weight = 1,
label = ~as.character(Lcontent),
clusterOptions = markerClusterOptions()) %>%
addPolygons(data = nbh_modified,
fillColor = "transparent",
color = "#000000",
fillOpacity = 0.8,
group = "Neighbourhood",
weight = 2) %>%
addPolygons(data = ctr,
fillColor = "transparent",
color = "#000000",
fillOpacity = 0.8,
group = "Census Tract",
weight = .8) %>%
addPolygons(data = disem,
fillColor = "transparent",
color = "#000000",
fillOpacity = 0.8,
group = "Disemination Area",
weight = .4) %>%
addLayersControl(overlayGroups =c("Neighbourhood", "Census Tract", "Disemination Area"),
options = layersControlOptions(collapsed=FALSE))